Set Up¶
# Import Libraries
import json # Reading and writing JSON
from glob import glob # File pattern matching
import os # Operating system access
import pathlib # File navigation
import holoviews as hv # Interactable plot
import panel as pn # Plot formatting
pn.extension() # Plot formatting
hv.extension('bokeh') # Interactable plot
import earthpy # Work with geospatial data
import xarray as xr # Multi-dimensional arrays
import hvplot.xarray # HVPlot and xarray
import rioxarray as rxr # Raster for xarray
import matplotlib.pyplot as plt # Plotting
import pandas as pd # Work with tabular data
import geopandas as gpd # Wowk with geospatial data
import hvplot.pandas # Plot geospatial data
Download Data¶
The dataset is saved as 'Gila River Vegetation.'
project = earthpy.Project(
'Gila River Vegetation', dirname='gric-vegetation-project')
project.get_data()
The 2020 American Indian Tribal Subdivision National (AITSN) borders from the 2020 US census are included in the data downloaded. More information about the shapefiles can be found here.
# Load in the boundary data
aitsn_gdf = gpd.read_file(project.project_dir / 'tl_2020_us_aitsn')
# Check that it worked
print(type(aitsn_gdf))
aitsn_gdf.head()
<class 'geopandas.geodataframe.GeoDataFrame'>
| AIANNHCE | TRSUBCE | TRSUBNS | GEOID | NAME | NAMELSAD | LSAD | CLASSFP | MTFCC | FUNCSTAT | ALAND | AWATER | INTPTLAT | INTPTLON | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2430 | 653 | 02419073 | 2430653 | Red Valley | Red Valley Chapter | T2 | D7 | G2300 | A | 922036695 | 195247 | +36.6294607 | -109.0550394 | POLYGON ((-109.2827 36.64644, -109.28181 36.65... |
| 1 | 2430 | 665 | 02419077 | 2430665 | Rock Point | Rock Point Chapter | T2 | D7 | G2300 | A | 720360268 | 88806 | +36.6598701 | -109.6166836 | POLYGON ((-109.85922 36.49859, -109.85521 36.5... |
| 2 | 2430 | 675 | 02419081 | 2430675 | Rough Rock | Rough Rock Chapter | T2 | D7 | G2300 | A | 364475668 | 216144 | +36.3976971 | -109.7695183 | POLYGON ((-109.93053 36.40672, -109.92923 36.4... |
| 3 | 2430 | 325 | 02418975 | 2430325 | Indian Wells | Indian Wells Chapter | T2 | D7 | G2300 | A | 717835323 | 133795 | +35.3248534 | -110.0855000 | POLYGON ((-110.24222 35.36327, -110.24215 35.3... |
| 4 | 2430 | 355 | 02418983 | 2430355 | Kayenta | Kayenta Chapter | T2 | D7 | G2300 | A | 1419241065 | 1982848 | +36.6884391 | -110.3045616 | POLYGON ((-110.56817 36.73489, -110.56603 36.7... |
Through some trial and error, I found latitude and longitude boundaries to focus the plot on the Gila River Indian Community. The AIANNHCE value 1310 can then be used to create a new GDF containing just the GRIC Boundaries.
# Generate plot of surrounding area with Tribal locations highlighted
aitsn_gdf.hvplot(
geo=True, tiles='EsriImagery',
frame_width=500,
legend=False, fill_color=None, edge_color='white',
hover_cols='all',
xlim=(-112.4, -111.4),
ylim=(32.8, 33.6)
)
WARNING:param.main: edge_color option not found for polygons plot with bokeh; similar options include: ['muted_color', 'bgcolor', 'line_color']
# Select and merge the subdivisions we want (1310)
gric_gdf = aitsn_gdf.loc[aitsn_gdf.AIANNHCE=='1310'].dissolve()
# Plot the results with web tile images
gric_gdf.hvplot()
# Create the boundary GDF
boundary_gdf = aitsn_gdf.loc[aitsn_gdf.AIANNHCE=='1310'].dissolve()
# Plot the results with web tile images
boundary_gdf.hvplot(
geo=True, tiles='EsriImagery',
fill_color='green', line_color='black',
fill_alpha=.3, # must be b/w 0-1
title='Gila River Indian Community',
frame_width=500)
Set up NDVI data¶
Create and combine a list of the raster files with the unique date for each raster now included as an extra dimension to the (x,y) spatial data.
# Get a sorted list of NDVI tif file paths
ndvi_paths = sorted(list(project.project_dir.rglob('*NDVI*.tif')))
# Display the first and last three files paths to check the pattern
ndvi_paths[:3], ndvi_paths[-3:]
([PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001145000000_aid0001.tif'),
PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001161000000_aid0001.tif'),
PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001177000000_aid0001.tif')],
[PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022209000000_aid0001.tif'),
PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022225000000_aid0001.tif'),
PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022241000000_aid0001.tif')])
# Loop through each NDVI image to add date dimension
doy_start = -25
doy_end = -19
ndvi_das = []
for ndvi_path in ndvi_paths:
# Get date from file name
doy = ndvi_path.name[doy_start:doy_end]
date = pd.to_datetime(doy, format='%Y%j')
# Open dataset
da = rxr.open_rasterio(ndvi_path, mask_and_scale=True).squeeze()
# Add date dimension and clean up metadata
da = da.assign_coords({'date': date})
da = da.expand_dims({'date': 1})
da.name = 'NDVI'
# Prepare for concatenation
ndvi_das.append(da)
len(ndvi_das)
154
Combine the list of rasters into a new data xarray "sliced" by the dates.
# Combine NDVI images from all dates
ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
ndvi_da
/tmp/ipykernel_1220/646390900.py:2: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly. ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date']) /tmp/ipykernel_1220/646390900.py:2: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly. ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
<xarray.Dataset> Size: 48MB
Dimensions: (date: 154, y: 203, x: 382)
Coordinates:
band int64 8B 1
* x (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
* y (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
spatial_ref int64 8B 0
* date (date) datetime64[ns] 1kB 2001-01-14 2001-01-16 ... 2022-01-24
Data variables:
NDVI (date, y, x) float32 48MB 0.8282 0.6146 ... 0.2146 0.2085# Plot first and last files side by side to visualize progress
first_ndvi = rxr.open_rasterio(ndvi_paths[0], mask_and_scale=True).squeeze()
last_ndvi = rxr.open_rasterio(ndvi_paths[-1], mask_and_scale=True).squeeze()
fig, axes = plt.subplots(1, 2, figsize = (12, 4))
first_ndvi.plot(ax=axes[0], cmap=plt.cm.PiYG, vmin=-1, vmax=1)
axes[0].set_title("Gila River NDVI - 2001")
boundary_gdf.plot(ax=axes[0], edgecolor='black', facecolor='none',
linewidth=1)
last_ndvi.plot(ax=axes[1], cmap=plt.cm.PiYG, vmin=-1, vmax=1)
axes[1].set_title("Gila River NDVI - 2022")
boundary_gdf.plot(ax=axes[1], edgecolor='black', facecolor='none',
linewidth=1)
<Axes: title={'center': 'Gila River NDVI - 2022'}, xlabel='x', ylabel='y'>
Calculate mean values for NDVI from the beginning of the dataset in 2001 to the water rights settlement in 2004 and an xarray of rasters for each year post 2004.
# Compute the mean NDVI values before water rights were restored
ndvi_pre = (ndvi_da
.sel(date=slice('2001', '2004'))
.mean('date')
.NDVI
)
ndvi_pre.plot()
<matplotlib.collections.QuadMesh at 0x7c4d6b87b850>
ndvi_post = ndvi_da.sel(date=slice('2005', '2022')).NDVI
ndvi_post_yearly = ndvi_post.groupby('date.year').mean('date')
print (ndvi_post_yearly)
len(ndvi_post_yearly)
<xarray.DataArray 'NDVI' (year: 18, y: 203, x: 382)> Size: 6MB
array([[[0.5380143 , 0.47134283, 0.4491 , ..., 0.17115715,
0.17115715, 0.175 ],
[0.53234285, 0.43379998, 0.44574285, ..., 0.14785714,
0.14785714, 0.16301428],
[0.50931424, 0.4504714 , 0.40725714, ..., 0.1275 ,
0.1056857 , 0.07895714],
...,
[0.18635714, 0.18635714, 0.17112856, ..., 0.17614286,
0.17614286, 0.17889999],
[0.17275715, 0.17275715, 0.16979997, ..., 0.21141426,
0.21141426, 0.2217714 ],
[0.18247142, 0.17312858, 0.18251428, ..., 0.20607142,
0.26507142, 0.2588 ]],
[[0.6633143 , 0.5207143 , 0.4350571 , ..., 0.14278571,
0.14278571, 0.14188571],
[0.6846143 , 0.5030286 , 0.35744286, ..., 0.11975714,
0.11975714, 0.12397142],
[0.5164143 , 0.42664286, 0.3745286 , ..., 0.13028571,
0.09927142, 0.07134286],
...
[0.14237143, 0.14237143, 0.13668571, ..., 0.16885713,
0.16885713, 0.17524286],
[0.1333857 , 0.1333857 , 0.13745713, ..., 0.21079998,
0.21079998, 0.21817143],
[0.15061429, 0.13785715, 0.13837144, ..., 0.16395713,
0.21935715, 0.2031714 ]],
[[0.6115714 , 0.45694283, 0.40044284, ..., 0.16761427,
0.16761427, 0.16799998],
[0.5378714 , 0.44529995, 0.35509998, ..., 0.15104285,
0.15104285, 0.15302856],
[0.33719996, 0.37427142, 0.4122714 , ..., 0.1349 ,
0.11072857, 0.09608571],
...,
[0.12699999, 0.12699999, 0.12062857, ..., 0.15535715,
0.15535715, 0.16347143],
[0.1187857 , 0.1187857 , 0.11717142, ..., 0.19084285,
0.19084285, 0.22217143],
[0.13205715, 0.11748571, 0.12167142, ..., 0.17574285,
0.19661427, 0.21658571]]], shape=(18, 203, 382), dtype=float32)
Coordinates:
band int64 8B 1
* x (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
* y (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
spatial_ref int64 8B 0
* year (year) int64 144B 2005 2006 2007 2008 ... 2019 2020 2021 2022
Attributes:
units: NDVI
AREA_OR_POINT: Area
18
# Calculate difference between mean ndvi before water rights and each year following the settlement
ndvi_diff = ndvi_post_yearly - ndvi_pre
# Check dimensions and coordinates for calculated difference
print(ndvi_diff.dims)
print(ndvi_diff.coords)
('year', 'y', 'x')
Coordinates:
band int64 8B 1
* x (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
* y (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
spatial_ref int64 8B 0
* year (year) int64 144B 2005 2006 2007 2008 ... 2019 2020 2021 2022
# Plot the difference
diff_plot = (
ndvi_diff.hvplot(x='x', y='y', cmap='PiYG', geo=True, groupby='year',
clim=(-0.5, 0.5),
title = 'Gila River Indian Community NDVI\n'
'Comparing 2001-2004 to Following Years')
*
boundary_gdf.hvplot(geo=True, fill_color=None, line_color='black')
)
diff_plot
BokehModel(combine_events=True, render_bundle={'docs_json': {'9151993b-ad93-44e6-af5e-51ce56d38647': {'version…
Move the slider to the bottom of the plot and set start to 2005.
years = list(diff_plot.kdims[0].values)
years.sort() # just to be safe
diff_plot_bottom = pn.panel(
diff_plot,
widget_location='bottom',
widgets={
'year': pn.widgets.DiscreteSlider(
name='year',
options=years,
value=years[0]
)
}
)
diff_plot_bottom
BokehModel(combine_events=True, render_bundle={'docs_json': {'2b12ca0f-b559-4be4-afb7-ff4da0790ccc': {'version…
# Save the new plot
diff_plot_bottom.save('ndvi_diff_plot_new.html', embed=True)
WARNING:bokeh.core.validation.check:W-1005 (FIXED_SIZING_MODE): 'fixed' sizing mode requires width and height to be set: figure(id='e0e0ac3a-719d-4a3e-aed1-9ed229c8e554', ...)